clear
est clear
cap mat define migrationfigure=J(105,10,.)
		local row=1
		

	use "$input/full_geoid_pair_panel.dta"
	
	

	replace migration=migration/(tot_pop/100000)

	//Make sure only include states that have a storm in at least one county at least once in our sample years
	bysort state_fips: egen max=max(storm)
		drop if max!=1
		
	//Create future pooled migration - total migrants in the subsequent five years
	sort geoid_pair move_year
	g post_storm=1 if storm==1|storm[_n-1]==1&geoid_pair==geoid_pair[_n-1]|storm[_n-2]==1&geoid_pair==geoid_pair[_n-2]|storm[_n-3]==1&geoid_pair==geoid_pair[_n-3]|storm[_n-4]==1&geoid_pair==geoid_pair[_n-4]|storm[_n-5]==1&geoid_pair==geoid_pair[_n-5]
	mvencode post_storm, mv(0)
	
	//Collapse to total migration from a county (not looking at paired county to county migration, just all migration out of a county)
	collapse (sum) migration (mean) share_approved totalapprovedihpamount,by(from_geoid move_year storm storm_fema post_storm state_fips)
	
	//Create IHS transformations
	ihstrans(migration)
	
	g total_storms=0
	//Create lags of storms (was there a storm in t-1 years)
	sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g storm_L`i'=storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
				replace total_storms=total_storms+storm[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
			
	label var post_storm "1(0-5 yrs post storm)"
	label var total_storms "All storms in t-5 years"
	label var storm "1(Storm year)"
	foreach var in 1 2 3 4 5{
		label var storm_L`var' "1(Storm year), t-`var'"
	}
	
	//Create leads of storms (was there a storm in t+1 years)
	sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g storm_E`i'=storm[_n+`i'] if from_geoid==from_geoid[_n+`i'] 
				 
		}	
			
	
	foreach var in 1 2 3 4 5{
		label var storm_E`var' "1(Storm year), t+`var'"
	}
	
	
	
	
	replace totalap=totalap/1000000
	
	forvalues var=1996(1)2013{
		g sample_`var'=(storm_L1==0&storm_L2==0&storm_L3==0&storm_L4==0&storm_L5==0&storm_E1==0&storm_E2==0&storm_E3==0&storm_E4==0&storm_E5==0&from_geoid==from_geoid[_n-5]&from_geoid==from_geoid[_n+5]&move_year==`var')
		mvencode sample_`var', mv(0) o
	}
		forvalues var=1996(1)2013{
		bysort from_geoid: ereplace sample_`var'=max(sample_`var')
	}
	
	local row=1
	local stack=1
	forvalues var=1996(1)2013{
		
		preserve
			keep if sample_`var'==1
			local high=`var'+5
			local low=`var'-5
			drop if move_year>`high'
			drop if move_year<`low'
			
			sort from_geoid move_year
			
			replace storm=storm[_n-1] if storm[_n-1]==1 & from_geoid==from_geoid[_n-1]
			
			
			
			if _N>0{
			
				reghdfe ihs_migration i.storm, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
				
					cap mat def migrationfigure[`row',1]=_b[1.storm]
					cap mat def migrationfigure[`row',2]=_se[1.storm]
					cap mat def migrationfigure[`row',3]=`var'
					bysort from_geoid: egen treated=max(storm)
					bysort from_geoid: g not_treated=(treated!=1)
					bysort move_year: egen sum=total(treated)
					bysort move_year: egen sum_non=total(not_treated)
					
					cap drop post
					
					
					sum sum, d
						local t1=r(mean)
					sum sum_non, d
						local t1_non=r(mean)	
						
					if `t1'>10{
						g treat_time=`var'
						g ttt=move_year-treat_time
						g post=ttt>=0
						g postXtreat=post*treated
						g stack=`stack'
						g stack_treated=`t1'
						g stack_control=`t1_non'
						save "$input/temp_v`var'.dta", replace
						local stack=`stack'+1
					}	
					drop sum
					
					duplicates drop from_geoid, force
					egen total=total(treated)
					sum total, d
						local t=r(mean)
						local c=_N-`t'
						cap mat def migrationfigure[`row',4]=`t'
						cap mat def migrationfigure[`row',5]=`c'
					local row=`row'+1
				
					
			} 
			else{
				
			}
		restore	
	}	
	
	graph set window fontface Arial
	
	
	clear

	forvalues var=1996(1)2013{
		cap append using  "$input/temp_v`var'.dta"
	}
	
	tostring ttt, g(ttt2)
		encode ttt2, g(ttt3)
		tab ttt, g(ttt_x)
		
	drop if stack_treated<10	
		
	egen total_treated=total( stack_treated)
		egen total_control=total( stack_control)
		g weight=((stack_treated/total_treated)/(stack_control/total_control))
		replace weight=1 if treated==1
	
	

	

	
	reghdfe ihs_migration ttt_x1##treated ttt_x2##treated ttt_x3##treated ttt_x4##treated ttt_x6##treated ttt_x7##treated ttt_x8##treated ttt_x9##treated ttt_x10##treated ttt_x11##treated [aweight=weight], absorb(i.from_geoid#i.stack  i.move_year#i.stack i.state_fips#c.move_year#i.stack) vce(cluster i.state_fips#i.stack)
		mat result = (r(table)["b", "1.ttt_x1#1.treated".."1.ttt_x11#1.treated"] \ r(table)["se", "1.ttt_x1#1.treated".."1.ttt_x11#1.treated"])'
		
		preserve
			clear
				svmat result, n(col)
				
					drop if se==.
					g n=_n
						*g odd=mod(n,2)
					*keep if odd==1
					g n2=_n
						replace n2=n2+1 if n>4
						
					insobs 1
						replace b=0 if b==.
						replace n2=5 if n2==.
					sort n2	
					
					g ul=b+2.67*se
					g ll=b-2.67*se
					
		twoway connected b n2 || rcap ul ll n2, xline(6) yline(0) legend(off) ytitle("% change in out-migration") xtitle("Years to treatment") xlab(1 "-5" 6 "0" 11 "5")|| rcap ul ll n2 if n2==6, color(cranberry)
			graph export "$figures/si6.png", as(png) replace
		restore